Stochastic Population Dynamics in Populations of Western Terrestrial Garter Snake with Divergent Life-histories

Comparative evaluations of population dynamics in species with temporal and spatial variation in life history traits are rare because they require long-term demographic time series from multiple populations. We present such an analysis using demographic data collected during the interval 1978-1996 for six populations of western terrestrial garter snakes (Thamnophis elegans) from two evolutionarily divergent ecotypes. Three replicate populations from a slow-living ecotype, found in mountain meadows of northeastern California, were characterized by individuals that develop slowly, mature late, reproduce infrequently with small reproductive Page 2 – Miller et al. – Stochastic population models effort, and live longer than individuals of three populations of a fast-living ecotype found at lakeshore locales. We constructed matrix population models for each of the populations based on 8-13 years of data per population and analyzed both deterministic dynamics based on mean annual vital rates and stochastic dynamics incorporating annual variation in vital rates. (1) Contributions of highly variable vital rates to fitness (λs) were buffered against the negative effects of stochastic variation, and this relationship was consistent with differences between the between the M-slow and L-fast ecotypes. (2) Annual variation in the proportion of gravid females had the greatest negative effect among all vital rates on λs. The magnitude of variation in the proportion of gravid females and its effect on λs was greater in M-slow than L-fast populations. (3) Variation in the proportion of gravid females, in turn, depended on annual variation in prey availability, and its effect on λs was 4to 23times greater in M-slow than Lfast populations. In addition to differences in stochastic dynamics between ecotypes, we also found higher mean mortality rates across all age classes in the L-fast populations. Our results suggest that both deterministic and stochastic selective forces have affected the evolution of divergent life-history traits in the two ecotypes, which, in turn, affect population dynamics. Mslow populations have evolved life history traits that buffer fitness against direct effects of variation in reproduction and that spread life-time reproduction across a greater number of reproductive bouts. These results highlight the importance of long-term demographic and environmental monitoring and of incorporating temporal dynamics into empirical studies of life history evolution. Page 3 – Miller et al. – Stochastic population models INTRODUCTION Life history theory focuses on understanding patterns of life history differences among species and the underlying processes that generate these differences (Stearns 1992, Roff 1992). Empirical examinations of population dynamics in the field have therefore focused on determining the ecological underpinnings of evolved differences in life histories. Central to ecological inference in this area has been the tradition of comparative studies within and among species on differences in mean demographic parameters (e.g., Promislow and Harvey 1990, Martin 1995, Reznick et al. 1996). Such comparisons of average demographic rates have yielded important insights about the causes of life-history differences such as how predation, a “topdown” agent of selection, and variable resources, a “bottom-up” selective force, are related to evolution in life-history traits (Martin 1995, Bronikowski et al. 2002, Walsh and Reznick 2009). In contrast, empirical examinations of temporal variation in population dynamics are much rarer, owing to limited availability of long-term data sets amenable to testing theoretical predictions (Morris and Doak 2004, Nevoux et al. 2010). Studies that incorporate temporal variability are necessary to extend the bridge between theoretical predictions and empirical observations that has served as the basis for our current understanding of life history evolution. Evolutionary models suggest that temporal variation in population vital rates and dynamics and in ecological interactions may impose strong selection pressure simultaneously on life history traits and demographic rates and should help explain patterns of life-history evolution (Orzack and Tuljapurkar 1989, Tuljapurkar 1990a, Benton and Grant 1999, Roff 2002, Koons et al. 2008, Horvitz and Tuljapurkar 2008, Tuljapurkar et al. 2009, Coulson et al. 2010, Tuljapurkar 2010). Methods used to study empirical population dynamics have recently been expanded to incorporate a stochastic demographic perspective into analyses (Morris and Doak 2002, Lande et Page 4 – Miller et al. – Stochastic population models al. 2003). Thus, since earlier notable reports (Benton et al. 1995), the empirical literature that utilizes stochastic models has been increasing (Coulson et al. 2001, Ezard et al. 2008, Morris et al. 2008). In addition, applied studies have incorporated stochastic dynamics into population models to guide conservation actions (Fieberg and Ellner 2001, Morris and Doak 2002, Lande et al. 2003) and to predict how changes in environmental variability will affect populations in the future (Boyce et al. 2006, Morris et al. 2008, Jonzén 2010). A rich set of theoretical and methodological resources exists that can improve insights from the analysis of long-term demographic data sets regarding the importance of temporal variation in a third area, that of life history evolution (Tuljapurkar 1990b, Levin et al. 1996, Pfister 1998, Tuljapurkar et al. 2003, Morris and Doak 2004, Haridas and Tuljapurkar 2005, Doak et al. 2005, Tuljapurkar et al. 2009). Specific outcomes from empirical analyses of long-term temporal dynamics are expected to advance our understanding of evolution in variable habitats (e.g., Morris et al. 2011). Important outcomes include testing the prediction that life histories should be buffered against the negative effects of variability (Pfister 1998), measures of the selection gradients associated with changes in stochastic variation in individual life history traits (Lande 1982, Caswell 2001, Tuljapurkar 2009), and measures of the direct relationship between year-to-year environmental variation and the fitness costs resulting from such variation (Davison et al. 2010, Tuljapurkar 2010). Fitness in temporally variable environments for density-independent populations is measured by the geometric mean of growth rates during individual time steps; therefore a cost may be imposed because higher variability decreases the geometric mean (Lewontin and Cohen 1969). Macro-level analyses have supported the view that selection has worked to buffer the fitness cost of variability by selecting for lower sensitivity to life history traits that are most temporally variable (Pfister 1998, Sæther and Bakke 2000, Gaillard and Yoccoz 2003, Dalgeish Page 5 – Miller et al. – Stochastic population models et al. 2010). Whether buffering of life-histories can explain differences among populations within species, and whether buffering occurs at the level of individual vital rates remain open issues. Comparisons of stochastic demography are also needed to test predictions about the sensitivities of different vital rates to changes in temporal variation (e.g., Ezard et al. 2008). Furthermore, investigations that identify the causes of temporal variation in vital rates (e.g., Coulson et al. 2001), and the ecological context which leads to differences in stochastic dynamics are particularly important in understanding how stochastic processes influence life history evolution (Koons et al. 2009, Tuljapurkar et al. 2009, Jonzén et al. 2010, Nevoux et al. 2010, Tuljapurkar 2010). Here we test these predictions about the importance of temporal variation viz. life-history evolution among spatially-distinct populations of the western terrestrial garter snake (Thamnophis elegans). Study populations in the vicinity of Eagle Lake in northeastern California originated from one ancestral source population that became differentiated into two genetically diverged ecotypes that have evolved distinctive morphologies and life history strategies (Bronikowski and Arnold 1999, Bronikowski 2000, Manier et al. 2007, Sparkman et al. 2007, Robert and Bronikowski 2010). Populations in this system correspond to either a lakeshore (Lfast) ecotype in which individuals grow fast, mature early, breed frequently with large litter sizes, and die young; or a meadow (M-slow) ecotype where individuals grow slowly, mature late, breed less frequently with smaller litter sizes, and live twice as long as L-fast individuals. These distinctive life-history strategies are maintained despite low levels of gene flow between the two ecotypes (Bronikowski 2000, Manier and Arnold 2005), and common environment experiments support the view that early life-history differences among populations have a genetic basis (Bronikowski 2000). In this report, we compare stochastic population dynamics for Page 6 – Miller et al. – Stochastic population models M-slow and L-fast populations to evaluate support for a relationship between temporal variability in environmental resources and associated population dynamics, and whether the relationship correspond to life history differences among populations. METHODS Study populations and field methods We collected data on six focal populations of western terrestrial garter snakes in the vicinity of Eagle Lake in Lassen County, California, representing three L-fast and three M-slow populations. From 1978-1996, demographic data were collected for these populations, corresponding to the lakeshore (L1, L2, and L3) and meadow populations (M1, M2, and M3). There was variation in the specific study years and length of study among populations and we report results from analysis of annual demographic models constructed for: 1980-1987 for L1; 1979-1988, 1994, and 1995 for L2; 1979-1988, 1994, and 1995 for M1; 1978-1988, 1994, and 1995 for M2; and 1980-1988, 1994, and 1995 for M3. A third lakeshore population (L3) was studied intensively between 1979 and 1988. However, limited recaptures of marked animals precluded estimates necessary to generate models with temporal dynamics for this population. Because the data provide an additional lakeshore data point for deterministic dynamics and for estimates of variance in reproductive parameters we include the population in these analyses. Snakes were captured during systematic searches of study areas. Areas were repeatedly searched on an annual basis from June to August with additional searches in May and October in some of the years. All snakes were individually marked by scale clipping or with a passive integrated transponder (PIT) tag. Snout-vent length was measured and sex was determined for all individuals. During the summer months, adult females were palpated to determine whether they were gravid and to estimate litter size based on the number of embryos. Accuracy of this Page 7 – Miller et al. – Stochastic population models method was validated by a strong correlation between litter size estimated by palpation in the field and litter size recorded at birth for the same individuals after being brought into captivity (n = 78, r = 0.88). To estimate survival rates, capture records for each individual were compiled, including recaptures in subsequent years, and these were used to generate unique capture histories for individuals throughout the study period (Lebreton et al. 1992). Anurans comprised an annually variable portion of snake diet compared to other prey, and presence of breeding anurans was closely tied to environmental conditions (i.e., precipitation and temperature; Kephart 1982, Kephart & Arnold 1982). In a previous study of the L2 population, Kephart and Arnold (1982) found that anurans comprised the majority of snake diet in years when a low proportion of captured snakes had empty stomachs, demonstrating a strong correlation between anuran presence and general food availability. When available, diet analysis has shown that snakes from both ecotypes eat the tadpoles and metamorphs of Pseudacris regilla and Bufo boreas but anurans compose a greater proportion of the diet in meadow habitats where fish are largely unavailable (Kephart 1982). Lakeshore snakes feed on Bufo boreas when they are available and in years where precipitation is especially high they can comprise more than half their diet. However, the majority of L-fast snakes diet is fish (Rhynichthys osculus) and leeches (Erpobdella spp.), which are abundant in lake habitats and predicted to be relatively consistent in availability among years. Based on these studies we made a number of predictions about how food availability would affect demography. We predicted that the presence of breeding anurans would be a strong predictor of annual variability for both populations, but that the availability of alternative food resources would mean that the effect on annual variability in lakeshore populations would be smaller than in meadow populations. We classified years during the study into those when breeding anurans were present, and tadpoles and metamorphs were abundant at Page 8 – Miller et al. – Stochastic population models study sites and those when they were not. We used detailed field notes about summer conditions, notes on anuran presence, and stomach contents from snakes captured in the field to classify each of the years. Predictions were tested using methods described below. Overview of Analyses We constructed a set of matrix population models which we used to gain insights into how stochastic population dynamics differed between the two life history ecotypes. We relied on a combination of deterministic and stochastic techniques as well as prospective and retrospective analyses (Caswell 2001). Our approach to model-development, parameterization, and analysis involved the following steps: 1) We defined the life cycle graphs and the corresponding Lefkovitch matrices (Lefkovitch 1965) we used for modeling the population dynamics for each of the populations. 2) We generated estimates of 5 vital rates needed to populate the models: the proportion of adult females that were gravid (PG), litter size (LS), neonate survival (SN), juvenile survival (SJ), and adult survival (SA). 3) For each of the populations we conducted a standard prospective matrix model analysis, estimating both deterministic and stochastic growth rates, sensitivities, and elasticities. This allowed us to compare differences in the means and variances of vital rates among populations and the relative strength of selection associated with potential changes in each of the means and variances. 4) We determined whether differences in life-histories among populations reduced the effects of variation (i.e., buffered) in vital rates on mean fitness (λs). We examined the relationship between the coefficient of variation and the elasticity of vital rates for Page 9 – Miller et al. – Stochastic population models each population, testing for an inverse relationship both within populations and within vital rates. 5) We examined the retrospective contribution of temporal variation and covariation in vital rates to annual variation in λ by analyzing a set of life-table response experiments (LTRE; Levin et al. 1996, Caswell 2000, 2001, 2010, Davison et al. 2010). We decomposed the total fitness costs related to annual variation in λ into contributions of variance in each vital rate and covariance among vital rates using a random design LTRE. We then decomposed fitness costs directly related to annual variation in the presence of breeding anurans into contributions from each vital rate using a stochastic LTRE. Life Cycle and Matrix Model In developing the population model used in our analyses we placed an emphasis on accurate estimation of temporal variation. To accomplish this goal we constructed a parsimonious structure for the model that allowed us to focus on resolving among-year variation of the components. For both ecotypes, we recognized three stages (neonates, juveniles, and adults). We assumed that the probability of being gravid and litter size were the same for all adult females. In addition, we assumed that transitions among stages happened at the same age for all individuals within a population. Individuals were assumed to spend 1 year in the neonate stage in both ecotypes. Based on previous work we stipulated that individuals spent either one year or three years as juveniles, for the L-fast and M-slow populations, respectively (Bronikowski and Arnold 1999). Finally, we were unable to estimate survival from birth (a latefall event) to the beginning of the first summer ( ). Because we were unable to estimate differences in this parameter among years and populations we used a fixed value of 0.8 , picked Page 10 – Miller et al. – Stochastic population models to approximate the median SN rate for a partial year. Consequently, variation in this parameter was not included in models. We constructed our models for the female portion of the population based on a prebreeding census design. We generated life cycle graphs and the corresponding annual transition matrices (At) and a mean annual matrix (Ā) as a function of stage specific survival ( ) and fecundity (F). (Figure 1). The life cycle analyzed for L-fast populations was a three stage model, 1, 2, and ≥3-year old, with each of the stages respectively corresponding with neonate, juvenile, and adults. For M-slow populations the life-cycle consisted of five stages: 1-year olds defined as neonates, 2, 3 and 4-year olds as juveniles, and ≥5-year old as adults. For both ecotypes, only adults contributed to the neonate stage through fecundity, which was estimated as a function of component vital rates, F = PG*LS*0.5* . The equation for fecundity has the assumption of a 50:50 sex ratio at birth, based on our many years of observed litter sex ratio. We were unable to estimate annual variability in survival for the L3 population, and so we only estimated the mean matrix, Ā. Thus, analyses for this population were limited to those based on deterministic growth rates and deterministic perturbation analyses. Estimation of Vital Rates We generated mean annual estimates and estimates of annual variance for the set of vital rates used in our population models, PG, LS, SN, SJ, and SA. Estimating variance based on raw annual estimates of vital rates will overestimate process variation in the rate because of the added variation due to sampling error. To account for this we used standard variance decomposition procedures to separate the sampling from process components of variation (Morris and Doak 2002). This separation was done post-hoc for LS and PG and as part of the estimation process using a hierarchical modeling approach for each of the survival rates. Estimates of PG were Page 11 – Miller et al. – Stochastic population models generated based on the proportion of adult females caught during June-August that were pregnant. LS was estimated from field data of litter size based on palpating adult females for young. We used methodology outlined by White (2000) that accounts for unequal sampling variance among years to estimate process variation in both PG and LS. Our estimation procedure for survival rates drew on standard methodology for mark-recapture estimates of survival (Lebreton et al. 1992) embedded within a hierarchical model and fit using Markov Chain Monte Carlo methods (White et al. 2008, Lukacs et al. 2008). This allowed us to directly estimate a random effect for annual variation in survival while accounting for sampling variation. Detailed description of survival estimation methods are found in Appendix A. Basic Matrix Model Analysis We focused on the influence of component vital rates vk rather than individual matrix elements aij in estimation of matrix sensitivities and elasticities (Morris and Doak 2004). We estimated the deterministic growth rate (λ1) and deterministic sensitivities and elasticities of vital rates ( ) and ) based on Ā for each population (Caswell 2001). Sampling error was calculated as a function of the standard errors of component vital rates. Deterministic sensitivities and elasticities of vital rates, vk, were calculated directly based on their relationship to the sensitivities for individual matrix elements, aij: ∑ , , ∑ , . Sensitivities measure the effect of absolute changes in vk on λ1 while elasticities of vital rates measured the effect of proportional changes in vk on log λ1. We examined population dynamics based on a stochastic population model for L1, L2, M1, M2, and M3. We calculated the stochastic growth rate (λs) using both Tuljapurkar’s Page 12 – Miller et al. – Stochastic population models (1990b) small variance approximation and stochastic simulation by randomly drawing among annual matrices, At (Morris and Doak 2002). The first method is approximate while the second is unbiased as the number of simulations increases. Both methods incorporate within year covariation among vital rates but assume that serial correlation among years does not occur. Initial results showed that differences in estimates of λs from the two methods were negligible for all populations and therefore we used Tuljapurkar’s approximation and analyses based on that approximation for all further analyses. The approximation draws on sensitivities estimated from the deterministic analysis for the m vital rates, log log ∑ ∑ , . (1) This equation demonstrates the effects of variation on fitness, showing that variance in a vital rate (var[vk] = cov[vi, vj], i=j=k) and positive covariance among the vital rates decreases λs, while negative covariance increases λs . These effects result in the constraint λs ≤ λ1. Sensitivity and elasticity of λs to changes in vital rates must account for changes in mean values and variance of the vital rates (Eq.1; Tuljapurkar et al. 2003, Doak et al. 2005, Haridas and Tuljapurkar 2005). Haridas and Tuljapurkar (2005) demonstrate that the sensitivity and elasticity of λs to changes in the means for vital rates are proportional to and from the deterministic analyses. Therefore, we focused on the effects of changes in temporal variation in each of the vital rates on λs. The effect on fitness of changes in variation for vital rates is of particular interest for stochastic analyses and sensitivities can be used as a direct measure of changes in fitness associated with an absolute change in variation in a vital rate. The small variance formulation for λs (Eq. 1) has been used to derive equations for sensitivities and elasticities of λs as a function of deterministic sensitivities (Doak et al. 2005, Haridas and Page 13 – Miller et al. – Stochastic population models Tuljapurkar 2005). Following the formulation by Doak et al. (2005) we calculated sensitivities of λs to the standard deviation of annual variation in a vital rate, , as ∑ , . (2) The stochastic sensitivity is not only a function of σi for each vital rate but also proportional to the correlation, ρij, of each vital rate to each of the other vital rates. The elasticity to λs was calculated as ∑ , . (3) Eq. 2 and 3 measure the effect of absolute and relative contributions of increased on λs based on the derivative of Eq. 1. Effects on λs include those due to increases in variances in the vital rate. In addition, increased variance in a vital rate will increase covariance even when correlation remains constant and Eq. 2 and Eq. 3 account for this associated effect. Note that Eσ and Sσ are only positive in cases in which negative correlations have a large contribution to Eq. 2 and 3.